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Abstract. The filamcntation instability driven by two spatially uniform and counter- 
streaming beams of charged particles in plasmas is modelled by a particle-in-cell 
(PIC) simulation. Each beam consists of the electrons and positrons. The four 
species are equally dense and they have the same temperature. The one-dimensional 
simulation direction is orthogonal to the beam velocity vector. The magnetic field 
grows spontaneously and rearranges the particles in space, such that the distributions of 
the electrons of one beam and the positrons of the second beam match. The simulation 
demonstrates that as a result no electrostatic field is generated by the magnetic field 
through its magnetic pressure gradient prior to its saturation. This electrostatic field 
would be repulsive at the centres of the filaments and limit the maximum charge 
and current density. The filaments of electrons and positrons in this simulation reach 
higher charge and current densities than in one with no positrons. The oscillations of 
the magnetic field strength induced by the magnetically trapped particles result in an 
oscillatory magnetic pressure gradient force. The latter interplays with the statistical 
fluctuations in the particle density and it probably enforces a charge separation, by 
which electrostatic waves grow after the filamcntation instability has saturated. 
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1. Introduction 

If a plasma is initially free of any net current, but if it has a nonequilibrium particle 
velocity distribution, then it can support the growth of magnetic fields through the 
redistribution of currents in space. This has been demonstrated first by Weibel [1] 
for a plasma with a bi-Maxwellian electron velocity distribution. The currents are 
rearranged through the growing plasma waves into filaments, which are separated 
by electromagnetic fields [2j. This is also the case for the instability driven by 
counterstreaming beams of charged particles, which is commonly referred to as the 
beam- Weibel instability or the filamentation instability (FI) [31 HI [5], [6] . 

The FI can generate strong magnetic fields in astrophysical environments such 
as the leptonic pulsar winds [7]. The FIs are also important for the generation of 
cosmological magnetic fields [8] and for inertial confinement fusion [9l [lO] , where laser 
pulses accelerate electron beams to relativistic speeds. Previous simulation studies have 
revealed various aspects of the growth and saturation of the FI. The FI driven by 
two counter-propagating beams of electrons has been examined by using particle-in- 
cell (PIC) and Vlasov simulations. Such studies have been performed in one spatial 
dimension (ID) [HI [121 [ISl [H [I5], in two spatial dimensions (2D) [31 [161 [HI [H] 
and in 3D [19]. The effects of a guiding magnetic field on counter-streaming electron 
beams have also been examined [El [20]. The simulation studies in Refs. [131 121] have 
investigated the impact of the ions on the non-linear stage of the FI. 

The FIs are usually triggered by the electrons. The ion filamentation is slower 
and often coupled through electrostatic fields to the electron filamentation. The FI 
involving only electrons couples to electrostatic fields during its quasi-linear growth 
phase. We refer to this source mechanism of electrostatic waves as the quasi-linear 
electrostatic instability (QEI). The electrostatic field amplitude grows in response to 
the QEI at twice the exponential rate of the magnetic field amplitude [121 [131 [13] and 
it oscillates around an equilibrium value after the FI has saturated. This equilibrium 
amplitude is such, that it excerts a force on the electrons that equals that of the average 
magnetic pressure gradient force (MPGF) [15]. This was demonstrated for the case 
of two counter-streaming, equally dense electron beams, for which the growth rate of 
the purely transverse FI is highest relative to those of the competing electrostatic two- 
stream instability and of the partially electromagnetic mixed mode instability [22]. The 
omission of wavevectors aligned with the beam velocity vector is thus most realistic. 

The pulsar winds carry with them positrons [7]. The impact of the positrons on 
the initial growth phase of the filamentation instabilities is understood, in principle, 
by solving the linear dispersion relation, as it has been done for example in Refs. 
[231 1211 [25] for the FI and in Ref. [26] for the Weibel instability. The nonlinear evolution 
of instabilities driven by beams of the electrons and positrons have been modelled with 
PIC simulations in 2D [23 [21] and in 3D [29]. 

In this paper, we consider two counterstreaming plasma beams, each of which 
contains an identical number of positrons and electrons. The FI is modelled in a short 
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ID simulation box, in which only one filament pair develops, as it has been done in Refs. 
[121 [ISl [IS] for the electron beams. The non-linear saturation of the FI is not captured 
correctly, since we exclude the filament merging and the multi-dimensional structure of 
the filaments [3]. However, the filaments are not circular if the beams are warm and if 
they have the same density. We find spatial intervals of the filaments that are planar 
over several electron skin depths [T71 [13 [30] . A ID simulation can give insight into the 
dynamics of such planar boundaries. The ID geometry allows us to freeze the filament 
pair just after the initial saturation and we can analyse the filaments in an almost time- 
stationary form. We can isolate a single filament pair to better understand its dynamics 
by omitting its collective interactions with the neighboring filaments. The restriction to 
one spatial dimension furthermore permits us to use a good statistical representation of 
the plasma phase space distribution and we can reduce the simulation noise. Accurate 
measurements of the fields and of the phase space distribution are thus possible. We 
choose initial conditions, which are identical to those in Refs. [Til [T5]. except for the 
positronic beam component that we include here. This allows us to compare directly 
the electromagnetic interaction of two counterpropagating electron beams with that of 
two counterpropagating beams of the electrons and positrons. 

We summarize our key results. The symmetry between the electrons and the 
positrons suppresses the QEI during the quasi-linear growth phase of the FI. The 
magnetic trapping model [H], which does not consider an electric field, accurately 
describes the saturation magnetic field in the simulation. The electrostatic field driven 
by the QEI would repel the electrons at the centre of the filament and attract those 
farther away [15], thus limiting the charge density accumulation due to the FI. Its 
absence implies that the filament confinement is stronger, if the positrons are present; 
higher peak density values can be reached and the spacing between the filaments is 
larger. The magnetic fields should reach higher spatial gradients compared to a system 
of counter-propagating electron beams due to the stronger currents. We could separate 
the purely magnetic FI driven by counter-propagating electron-positron beams from a 
secondary electrostatic instability (SEI), which we show to be unrelated to the QEI. 
The waves the SEI drives have a broadband wavenumber spectrum. We bring forward 
evidence for a connection between the SEI and the spatio-temporal oscillations of the 
MPGF. These oscillations occur on time scales, which are comparable to the inverse 
plasma frequency, and on spatial scales of the order of a Debye length. The MPGF can 
thus interplay through these oscillations with the statistical fluctuations of the plasma, 
by which a charge separation can occur. We propose that this separation breaks the 
initial symmetry of electrons and positrons and destabilizes the filament. The power of 
the waves driven by the SEI, which have a frequency that is close to the plasma frequency, 
grows in time. Its power can be fitted as a function of time by two exponential functions, 
which are separated by a break in the growth rate. 

The manuscript is organized as follows. Section 2 briefly describes the PIC 
simulation method and it discusses our initial conditions. Section 3 presents our 
simulation results, which are then discussed in Section 4. 
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2. The PIC simulation method and the initial conditions 

The standard particle-in-cell (PIC) method [31] can model the processes in a coUisionless 
kinetic plasma. It approximates the plasma phase space distribution by an ensemble of 
volume elements or computational particles (CPs). The ensemble properties of the CPs 
are an approximation to the ensemble properties of the corresponding physical plasma 
species. Each CP with index i of the species j can have a charge qj and mass mj that 
differ from those of the plasma species they represent, e.g. the mass rUe and the charge 
— e of an electron. The charge-to-mass ratio must be preserved though. 

The CPs follow trajectories in the simulation domain, which are determined by the 
electric E and the magnetic B fields. The electromagnetic fields and the global current 
J are defined on a spatial grid. The equations that are solved by the PIC code are 



The electromagnetic fields are evolved in time using the Faraday law and the Ampere 
law (Eq. ([1])). Equations ([2]) are typically fuUfiUed as constraints, or they are enforced 
by correction steps. The virtual particle method [32] our code is using is fullfiUing 
Poisson's equation as a constraint and V ■ B = is solved exactly in ID and to round- 
off precision in higher dimensions. The trajectories of the CPs are updated using the 
Eqns. ([3]). We refer the interested reader to the Refs. [3ll [321 [33l [SI [351 EB] for a more 
thorough discussion of numerical PIC schemes. 

Equations (l)-(3) can be scaled to physical units with the total plasma frequency 
i^p = {nte^ /meeo)^^"^ and the skin depth Ae = c/up, where the total particle number 
density Ut = J2j nj is summed over the four leptonic species, which are equally dense. 
The quantities in physical units denoted by the subscript p are obtained from the 
normalized ones by substituting Ep = tUpCmgE/e, Bp = Upm^'B/e, Jp = ecutJ, 
Pp = eutp, Xp = AgX, tp = t/up, Vp = vc and pp = mgCp. The charge qj, in this 
normalization, is 1 for the positrons and -1 for the electrons, while ruj = 1. We also 
normalize Q = u/ujp and k = kpc/up, where u, kp have physical units. 

Both beams in our simulation study move in opposite z-directions with the speed 
modulus Vb = 0.3, giving the relative beam speed 2vb/ (1 + v"^) ^ 0.55. This Vb is about 
half that used in Ref. [28]. Each of the two beams consists of one electron species 
and of one positron species and all four species are equally dense. Initially all beams 
are spatially uniform. The velocity distribution in the rest frame of each beam is a 

1 /2 

Maxwellian with the thermal speed Vt = [kT/nif,) ' = Vb/lS in all directions. The ID 
simulation box with its periodic boundary conditions is aligned with the x-direction. 
Only waves with wavevectors parallel to x can grow and we use the scalar wavenumber 
k. The length L = 1.25Ae of the simulation box is identical to that of the shortest one in 
Ref. [15], if we neglect the positron contribution to Ag. The simulation box is subdivided 



VxE=--. VxB=- + J, 
V ■ B = 0, V ■ E = p, 

— ^ = qj (E[xi] + Vi X B[xi]) , Pi = nij Vj r(vi) 




(1) 
(2) 

(3) 



Numerical studies of the electron-positron filamentation instability in ID 



5 



M 10' 



1 10- 

LU 



CD 
T3 



>, 



10 



10' 



-4 



,-2 



,-6 



-8 








50 



100 



150 



Time t 



Figure 1. The energy densities in units of the total energy in the simulation: The 
uppermost curve Eby corresponds to By. Its exponential rate, which is twice the fli 
of By, is w 0.5 during 10 < i < 50. It oscillates around an equilibrium after that. The 
curves of Eez and Eby grow initially at the same exponential rate, but Eez decreases 
for t > 50. The electrostatic Eex starts to grow at i « 45, when Eby has saturated. 



into 500 grid cells of equal length A^.. The ratio A^^Up/vt = 0.15 and the Debye length 
is resolved well. Each of the four plasma species is represented by 2.45 x 10^ CPs. The 
E = and B = at the simulation's start. The total simulation time is tsim = 177, 
which is subdivided into 10^ time steps A^. 

3. Simulation results 

The FI driven by a plasma flow along z and a simulation box that is aligned with x will 
lead to the initial growth of a magnetic field along y. The growing net current Jz{x) 
will also result in a growing by Ampere's law. An electrostatic field would grow 
in the case of a system that is composed of two electron beams. 

Figure [1] displays the energy densities of these field components, which we denote as 
Eby, Eez and Eex- The energy density Eby of the magnetic By component dominates 
and it reaches a few percent of the total energy, in line with previous simulations |28j . 
The exponential growth rate of By is Qi ~ 0.25, which is close to the expected value 

1 /2 

~ ffe/7f, for cold beams. The Eez grows at the same rate as Eby, but its values 
are three orders of magnitude less. The energy density Eex of the electrostatic field 
grows after the FI has saturated. The growth rate of Eex oc E"^ can be fitted with an 
exponential function with the growth rate ^ 0.13 between 45 < t < 90 and with 
a second, slower growing one for t > 100. The growth of Eex in Fig. [1] can not be 
attributed to the QEI, because then the Eex should grow at twice the exponential rate 
of Eby until t ~ 50 and oscillate around an equilibrium value after that time. 

We look in more detail at the fields to better understand the cause of the oscillations 
of Eby around its equilibrium and the source mechanism of the growth of Eex- The 
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Figure 2. (Colour online) The electromagnetic fields: Panels (a) and (b) show By 
and Ez, respectively. Both are correlated and their phase difference is 90°. The Ez 
field is damped in time. Panel (c) shows the spatial power spectrum logioP{k, t) of By. 
The ki mode dominates, but harmonics can be seen. Panel (d) shows the frequency 
power spectrum logiQP{ki, f2) of By. Peaks are found at il = and fl « 0.3. Panel (e) 
evidences a logiQP{k,t) of Ez resembling that of By in (c) and the logj^g P(fci, f2) of 
Ez in (f) also has a maximum at 17 ?» 0.3. The power in (c,e) and (d,f) is normalized 
to the same value. 



field components driven by the FI are investigated in Fig. [2l The By, saturate at 
t ~ 50. The By{x,t) then remains practically stationary, while Ez{x,t) is damped. 

We compute the spatial power spectrum P{k,t) of By,Ez and E^ by a ID Fourier 
transform over the full box length L and by taking the square of the amplitude modulus. 
The frequency power spectra (dispersion relation) are obtained through a 2D Fourier 
transform of the field data over the box length L and for t > 50. The amplitude moduli 
are then squared to give P{k,fl). The base-10 logarithm of P{k,t) of By in Fig. EJj^c) 
evidences, that most power is concentrated in the mode ki, with kg = 27is/L, and that 
the power in this mode is oscillating. Weaker harmonics with uneven s also occur. The 
P{ki,Q) of By reveals peaks aX Q = and Q ~ 0.3. The peak at Q = dominates, 
because By is practically stationary after t = 50. The P{k, t) and the P{ki, Q) of E^ in 
the Figs. [2t^e,f) resemble qualitatively those of By, but they are weaker and the peak 
with f2 = is absent in the i^^-field. 

The Ex{x,t) in Fig. [3] reveals no spatial correlation with By{x,t). The P{k,t) of 
the E^-field in Fig. [3](b) and of the i^j^-field in Fig. [2]^c) show no link between the k of 
the dominant waves. The spatial power spectrum P{k, t) in Fig. [3]^b) evidences instead 
a wave growth over a wide band of k. The E^ can not be driven by the QEI, because 
then the MPGF, which is OxPby = BydBy/dx in our normalization, would imply that 
Ex = Q whenever By = Q oi dBy/dx = and that it should oscillate in the /c2 mode. 
The P{k, Q) of E^ in Fig. [3](c) reveals that the strongest waves have a f2 ^ 1. 

The particle phase space distributions provide more information. The electrons of 
the beam 1 (moves in positive z-direction) are species 1 and the positrons are species 3. 
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Figure 3. (Colour online) The electrostatic field: Panel (a) shows Ex{x,t). 
Fluctuations grow at late times, which are distributed seemingly randomly in space. Its 
logj^Q P{k, t) in (b) evidences a broadband spectrum, that shows no dominant k. The 
logj^Q P{k, f2) peaks at f2 w 1 and A: > 0. The strongest fluctuations are propagating 
waves. The P{k,t) and P(fc,r2) are normalized to the corresponding values of By. 




(d) Position X (e) Position X (fj Position X 



Figure 4. (Colour online) The base-10 logarithm of the phase space densities at the 
time t = 50 in units of CPs: Panels (a-c) display the f{x,Vz) of the species 1, 2 and 
4 (positrons), (d-f) show the f{x,Vx) corresponding to the panels above them. The 
distributions in (d,f) are practically identical and (a,c) can be mapped onto each other 
by changing the sign of Vz. The f{x,Vx) of species 1 and 2 are shifted by L/2. 



The electrons and positrons of beam 2 are denoted as species 2 and 4, respectively. The 
current fluctuations of PIC simulations imply that Bz,Ey ^ 0. These fields correspond 
to waves that propagate in form of the high-frequency electromagnetic modes. However, 
the peak value of Eez exceeds the energy densities of Ey, by 6 orders of magnitude 
and the latter remain practically constant during the simulation time. We can restrict 
our investigation to the phase space projections f{x,Vx) and f{x,Vz), because = 0, 
Bz ^ and Ey ^ 0. The f{x,Vz) will reveal the electromagnetic structures, while the 
electrostatic structures are represented by f{x,Vx)- Figure H] shows the phase space 
distributions of the species 1,2 and 4 at t = 50. The electrons of beam 1 and 2 are 
separated in space, which is typical of the FI driven by counter-streaming electron 
beams. The species 4 shows an f{x, Vx) that is identical to that of species 1. The f{x, Vz) 
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0.2 0.4 0.6 0.8 1 1.2 

(b) Position X 

Figure 5. Panel (a) The number density distribution Nj{x,t = 50) of the four 
species j in units of their mean density A^o- The curves Ni{x) and N4{x) match 
(their differences amount to less than the thickness of the curves) and they peak at 
X = 0. The also matching N2{x) and N'i{x) have their maximum at x = L/2. Panel 
(b) plots AA^ = qjNj{x) / ANq, with qj = 1 and —1 being the positron and electron 
charge. The fluctuations amount to less than fO^^ and they are by their scale size 
K, A-c the statistical fluctuations due to the finite number of CPs per cell. 

distribution of species 1 can be mapped onto that one of species 4 by switching the sign 
of Vz- Both these observations are expected, because in the absence of a significant 
the opposite speed and charge of both species cancel. The Lorentz force thus displaces 
both species in the same way. The relation is also the same between the species 2 and 
3. The high degree of symmetry is also demonstrated by the movie 1, which animates 
in time the projected phase space distributions of species 1, which is added to that of 
species 4. The color scale shows the base-10 logarithm of the number of CPs. 

Figure [5] displays the number densities Ni{x) = J fi{x, Vx)dvx for each of the species 
i, normalized to the mean value Nq = {Ni{x)),^. The A''i = and = ^3 within 
the resolution of the image. The charge density modulations reach a peak value ~ 4A'o. 
This peak value is higher by a factor three and, accordingly, the filament confinement is 
thus stronger here than in Ref. [15], which did not consider positrons. Each species is 
represented by iVp = 4.9 x 10^ CPs per cell. The statistical fluctuations of the particle 
number are thus Np~^'^ ^ 5 x 10~^, which is comparable to the measured charge density 
fluctuations AA^ = {4:Nq)~'^ qiNi < 10~^ in Fig. Mjo)- This observation, together 
with that the AA^ oscillates on a scale Ax, implies that the fluctuations at t = 50 
are due to the finite number of CPs per cell. These fluctuations have a fc-spectrum 
that is qualitatively similar to that of thermal noise. The fluctuation amplitude AA^ is 
increased within the filaments, e.g. in the interval 0.5 < x < 0.75, because we have not 
normalized it to the local particle number, but to the average density. 

Figure E] displays the phase space distributions of the species 1, 2 and 4 at t = 177. 
The f{x,Vz) are qualitatively unchanged compared to those in Fig. [Hand species 1 and 
4 are still symmetric to a change of the sign of Vz- The distributions f{x,Vx) reveal 
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Figure 6. (Colour online) The base-10 logarithm of the phase space densities at the 
time t = 177 in units of CPs: Panels (a-c) show f{x,Vz) of the species 1, 2 and 4, 
respectively, (d-f) show the f{x,Vx) corresponding to the panels above. Species 1 and 
4 reveal similar f{x,Vx) and distributions f{x,Vz) that are qualitatively symmetric to 
a change of the sign of Vz- The /(x, Vx) of species 1 and 2 are shifted by L/2. 
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Figure 7. (Colour online) Panel (a) The Nj{x) at the time t=177 of the four species 
j in units of A^o- The electron distributions Ni{x) and N2{x) are denoted by the 
solid curves and N3,N4 are dashed and blue. The curves Ni{x) and N4{x) peak at 
X = 0. The curves for electrons and positrons do not match. Panel (b) plots the 
AA^ = J2j Qj^ji^) /4:No, with qj = 1 and —1 being the positron and electron charge. 
The fluctuation amplitude is about 10 percent and they oscillate on scales Ax. 



that the particles are concentrated at the same positions as in Fig. HI The particles 
have been heated up along the x-direction in between the dense filaments. The small 
scale structures, in particular in f{x,Vx), differ for species 1 and 4 and we expect now 
clear charge density modulations. Figure [7] compares Ni{x) with N4^{x) and N2{x) with 
N^lx). Their differences can now be seen even from the Ni. The value A^4(x = 0.08) 
exceeds that of A'^i(x = 0.08) by about 0.4 and fluctuations of this size occur in the 
entire box, which is demonstrated by the AA^ in Fig. [7](b). The typical length scale of 
the oscillations is about 0.1 3> Ax and they are caused by the waves in Fig. [3l 

The peak amplitude of Er^ in Fig. [3] is ~ 5 x 10"'^ at t = 177. The Lorentz 
force VbBy for Vb = 0.3 and a maximum By f» 0.07 at t = 177 (Fig. [2]) is larger than 
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Figure 8. (Colour online) Panel (a) shows the magnetic pressure gradient force 
BydBy/dx, which oscillates twice as fast in space as By. The MPGF is not stationary in 
time. Panel (b) displays the base-10 logarithm of the modulus of the spatial frequency 
spectrum A{x, il) of the MPGF. It is normalized to its peak value at = 0. Strong 
oscillations close to a; = i/4 and x = 3L/4 reach = 1. 

the electrostatic force by a factor of four. The electrostatic force is weaker before this 
time, while By is constant after t = 50. Thus, the electrons experience for most of the 
simulation duration the only as a perturbation and the gyrofrequency of the CPs is 
determined by By and Vz ~ Vb, which is also demonstrated by the movie 1. The rotation 
of the dense filament in the distribution f{x, v^) has an approximately constant angular 
velocity in phase space. The structures in f{x,Vz) and, thus, will change with this 
characteristic frequency and impose the modulation of By and with fl ?s 0.3 in the 
Fig. [2]|^d,f). This frequency is ~ Qi, and consistent with the magnetic trapping. 

The charge density fluctuations in Fig. Mjo) at t = 50 were at noise levels, which 
was expected from the low Eex in Fig- [U at this time and no QEI occured. It is only 
after the FI has saturated, that the electrostatic waves grow due to the SEI. The SET 
yields the growth of Eex in the interval 50 < t < 177 and the waves in Fig. [3] have no 
obvious correlation with the By in Fig. [21 We will now identify one potential cause of 
the SEI. The force imposed by a magnetic pressure gradient on a current J is 

J X B = -VBV2. (4) 

Only By is growing in our ID simulation box and the only possible spatial derivative 
is dx- Equation m simplifies to JzBy = BydBy/dx. The By is strong for t > 50 and 
dBy/dx 7^ 0. The MPGF on the right-hand side does, therefore, not vanish and the 
particles are accelerated. 

Figure [8] shows the MPGF after t = 50 and A{x,Q), which is its spatial frequency 
spectrum obtained by a Fourier transform over time. The maxima and minima of dxPsY 
are stationary in space but its magnitude oscillates in time. The time average of the 
MPGF is positive to the right of the positions a; = and x = L/2 and negative to the 
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left of these positions. The x = and x = L/2 coincide with the stable equilibrium 
points of the respective filaments. The MPGF is thus accelerating the particles away 
from the equilibrium points. The MPGF is, however, weaker than the drift force qjVbBy, 
which is responsible for the filament confinement (magnetic trapping). 

Consider the filament formed by the species 2 and 3 at t = 50, which is centred 
at X = L/2 in Fig. [5] and Fig. ^a). The electrons of species 2 move with the velocity 
~ —VbZ, while the positrons of species 3 move with ^ vi,z. Their currents have the 
same sign and the MPGF accelerates the electrons and the positrons into the same x- 
direction. As long as the positrons and the electrons have the same density everywhere, 
the Ex does not grow, because J^^ = in dtE^ + Jx = 0. The term dyB^ — dzBy = 
in the ID geometry. The MPGF does thus not result here in the Ex-G.e\d discussed in 
Ref. [T3l [15], which has a wavelength that is half of that of the wave in By. However, 
the finite number of CPs introduces statistical fluctuations in the charge density (Fig. 
|5]) that imply that the MPGF accelerates locally (on Debye length scales) a different 
number of electrons and positrons, by which 7^ 0. An electric field grows, that tries 
to restore the charge neutrality. The MPGF can couple to these fluctuations, if the 
force gradient is high (comparable to the spatial scale of the fluctuations) and if the 
force oscillates with the characteristic frequency of the fluctuations. Their dominant 
oscillation frequency in the spatial intervals with ^ is the plasma frequency and it 
is the upper- hybrid frequency otherwise. Figure [^c) shows that the growing fluctuations 
have a broad frequency band and peak at ~ 1. The spatial gradients of the plasma 
density will, however, influence the fluctuation spectrum [37] . A force interacting with 
such fluctuations should thus also have a broad frequency band. 

The spectrum A{x, Q) in Fig. Mjo) reveals that the highest frequencies can be 
reached by the MPGF at x ~ 0.1, 0.5, 0.7 and at 1.1, although not with a high driving 
amplitude. The positions close to x = L/4 and x = 3L/4 experience stronger oscillations 
up to = 1. The MPGF changes also on spatial scales comparable to the Debye length, 
which is Vt/c = 1/60 in our normalization. An example is here a change at t ~ 70 
and X ~ 0.325 in Fig. [8](a). The MPGF changes from 0.02 to -0.02 over a distance 
~ 0.03. The connection between the growth of Eex and the oscillations of the MPGF 
is evidenced also by a comparison of Fig. [T]with Fig. [H](a). The growth of Eex slows 
down at t ~ 100. The oscillations in Fig. [8](a) close to L/4 and 3L/4 are more intense 
and shorter in duration before t = 100 and the spatial gradients are higher. The MPGF 
can thus couple easier energy to the charge density fluctuations. The spatio-temporal 
oscillations of BxPey soften up after t = 100, and the growth of Eex slows down. 

4. Discussion 

In this paper we have considered the filamentation instability (FI) driven by two 
counterstreaming beams, each consisting of the electrons and positrons. The beams 
have initially been spatially uniform and the electromagnetic fields were set to zero. The 
ID simulation is aligned with the x-direction, which excludes the merging of filaments 
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beyond a certain size and limits the physical realism of its nonlinear evolution [3]. 
However, the filaments are usually not circular but elongated pn[T8| [30]. A ID geometry 
may thus be a valid approximation for those parts of filaments, which are quasi-planar 
such as the one investigated in Ref. [18]. The two beams move in opposite z-directions 
at the same speed modulus Vb = 0.3, which is sufficiently low to exclude significant 
relativistic effects. It is sufficiently high to obtain a linear growth rate of the FI that 
is comparable to those of the electrostatic modes and of the mixed modes. This may 
be further aided by the equal beam densities, which favor the filamentation instability 
if positrons are absent [22]. The strong electric fields with their oblique polarization 
in the PIC simulation in Ref. [28] demonstrate, however, that in particular the mixed 
modes would compete with the FI in a more realistic 2D or 3D simulation. 

Our aim has been to obtain further insight into the dynamics of a filament pair 
formed by two counterpropagating beams of the electrons and positrons and to compare 
it with that of an electron filament pair. The simulation parameters are identical to those 
in the simulation in Ref. [15] that did not take into account the positrons. The short 
simulation box allows us to use a good statistical plasma representation and low noise 
levels. The charge density fluctuations inherent to the PIC simulation method provide 
noise over the full band of wavenumbers resolved by the simulation, out of which the 
wavemodes and any secondary instabilities can grow. 

The FI redistributes the beams of charged particles in space into current filaments 
[3]. Our simulation box length L and the periodic boundary conditions allow only a 
single pair of filaments to grow. The centres of the electron filaments are spatially 
separated by the distance L/2 and this is also the case for both positron filaments. The 
filament formed by the electrons of one beam coincides with the filament containing the 
positrons of the second beam. Their phase space distributions f{x,Vx), which hold the 
information about the electrostatic structures, match and their phase space distributions 
f{x, Vz) can be mapped from one to the other by a change of the sign of Vz- The currents 
and densities of both components add up. The symmetry between the electrons and 
positrons within the same filament implies that the MPGF accelerates them into the 
same direction. We have found that their partial currents cancel each other until the FI 
saturates, implying that no can grow due to the QEI. 

The simulation in Ref. [I5] evidenced that this would accelerate the electrons 
away from the centre of the filament. The presence of the positrons in the simulation 
in the present work removes this repulsion. Higher charge and current densities can be 
reached by electron-positron beams compared to those containing only electrons and 
the spacing between the filaments is larger. The magnetic fields can reach higher spatial 
gradients. The magnetic fields grew to an amplitude set by magnetic trapping jllj . 

We have identified a SEI that is leading to a growing electrostatic energy density 
after the FI saturated, resulting in a broadband wave spectrum. The exponential growth 
rates of these waves remained well below that of the FI and they show no correlation 
with the MPGF. These observations are conflicting with the properties attributed to 
the QEI, implying that the SEI and the QEI must have different source mechanisms. 
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The finite number of computational particles imply statistical fluctuations of the 
charge density on Debye length scales. We have found here that the oscillation spectrum 
of the MPGF involves a wide band of frequencies that reach a maximum that exceeds 
the plasma frequency and that the MPGF changes significantly on scales comparable 
to the Debye length of the plasma after the FI has saturated. The MPGF can thus 
couple to the statistical fluctuations of the charge density. We have proposed that this 
coupling amplifies the fluctuations of and E^^ and that this is the cause of the SEI. 

This hypothesis can be tested and the properties of the SEI can be examined 
in more detail with Vlasov simulations. They solve the Vlasov-Maxwell equations 
directly and do not approximate the plasma by phase space blocks. They are thus 
free of noise due to statistical plasma density fluctuations. The SEI should not develop 
in Vlasov simulations, unless seed fluctuations are introduced. The filaments, which 
are plasma structures resulting out of magnetic instabihties, would be more stable in 
Vlasov simulations than in PIC simulations. Previously, such an enhanced stability has 
only been reported for nonlinear structures (electron phase space holes) evolving out of 
purely electrostatic Buneman instabilities [38] • A localized charge density perturbation 
can be introduced into a Vlasov simulation and its interplay with the MPGF can be 
investigated. Such studies can not be performed with PIC codes, where we always have 
an ensemble of charge density perturbations interacting with the MPGF. 

The energy density Eex of the waves driven by the SEI remained below 10~^ E by- 
The -Ea;-field will, however, not be negligible after t = 177, because the electric force qjE^ 
competes with the Lorentz force qjVzBy and Vz < 1. The ratio between the strongest 
Lorentz force and the strongest electric force has been 4 at t = 177 and the instability 
is growing. The SEI is also likely to play a more important role in PIC simulations 
that use only a low number of particles per cell. The higher relative charge density 
fluctuations imply, that the SEI grows from higher initial amplitudes, which reduces its 
growth time. The initial electric noise power scales approximately inversely proportional 
to the number of particles per cell. A reduction from our 4.9 x 10^ particles per cell to 
50 particles per cell will increase the initial power of E^ by the factor 10'^. 

The energy density of the electric field has grown to a significant fraction of that 
of the magnetic field in the ID and 2D simulations of counterpropagating electron 
beams [121 CSl [HI [iSl [E] , and it should influence the interplay of the filaments during 
their nonlinear evolution. No significant electric fields are observed during the quasi- 
linear growth phase of the FI and immediately after its saturation, if each beam carries 
electrons and positrons with an equal density. Our future studies will thus assess the 
impact of an absent electric field on the filament size distribution. This distribution can 
be approximated by a Gumbel distribution in ID if no positrons are present [13]. We 
will also examine with 2D simulations how the filament size distribution orthogonal to 
the beam velocity vector evolves in time when positrons are present. The characteristic 
filaments size increases in the absence of the positrons linearly with time [17]. The 
impact of the electrostatic and mixed modes on the evolution of the filaments will also 
be addressed by 2D simulation studies that contain the beam velocity vector. 
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